using System;
using System.Collections;
using System.Collections.Generic;
using System.Text;
using Data;
using Data.Probabilities;
using System.IO;
using System.Text.RegularExpressions;
namespace Logic.IB
{
    class DistanceFunction
    {
        ProbabilityData probData;
        private static double[,] distanceMatrix;
        private static bool isInitialized = false;
        private ArrayList aLocations;
        private ArrayList bLocations;

        public DistanceFunction(ProbabilityData probData)
        {
            this.probData = probData;
            aLocations = new ArrayList();
            bLocations = new ArrayList();
            distanceMatrix = new double[probData.NumberOfGenomes, probData.NumberOfGenomes];
            for (int i = 0; i < probData.NumberOfGenomes; i++)
                for (int j = i + 1; j < probData.NumberOfGenomes; j++)
                    distanceMatrix[i, j] = dist(i, j);
            isInitialized = true;
        }


        public DistanceFunction(StreamReader matrixFile)
        {
            int numberOfGenomes=0;
            Regex reg = new Regex("[0-9]+");
            MatchCollection m = reg.Matches(matrixFile.ReadLine());
            //Checks the number of genomes;
            if (m.Count != 1)
            {
                matrixFile.Close();
                throw new Exception("Bad Format in distance matrix file");
            }
            else
            {
                numberOfGenomes = Int32.Parse(m[0].Value);
                distanceMatrix = new double[numberOfGenomes, numberOfGenomes];
                for (int i = 1; i<numberOfGenomes; i++)
                {
                    reg = new Regex("([0-9])+[.][0-9]+");
                    m = reg.Matches(matrixFile.ReadLine());
                    if (m.Count != numberOfGenomes)
                    {
                        matrixFile.Close();
                        throw new Exception("Bad Format in distance matrix file");
                        System.Windows.Forms.MessageBox.Show("format error");
                    }
                    else
                    {
                        for (int j = 1; j < i; j++)
                        {
                            distanceMatrix[j, i] = Double.Parse(m[j].Value);
                        }
                    }
                }
                matrixFile.Close();
                isInitialized = true;
            }
        }
        

        public ArrayList getMatrixInArrayList()
        {
            if (isInitialized == true)
            {
                ArrayList matrix = new ArrayList();
                string line = "";
                for (int i = 0; i <probData.NumberOfGenomes; i++)
                {
                    line = "";
                    for (int j = 0; j < probData.NumberOfGenomes; j++)
                    {
                        if (j == i)
                            line = line + "   0.0";
                        else
                        {
                            if (j < i)
                                line = line + "   " + distanceMatrix[j, i];
                            else
                                line = line + "   " + distanceMatrix[i, j];
                        }
                    }
                    matrix.Add(line);
                }
                return matrix;
            }
            else
            {
                throw (new Exception("You need to initialize the matrix first"));
            }
        }

        public static double getDistance(int genomeA, int genomeB)
        {
            if (isInitialized)
            {
                if (genomeA == genomeB)
                    return 0.0;
                return (genomeA < genomeB) ? distanceMatrix[genomeA, genomeB] : distanceMatrix[genomeB, genomeA];
            }
            return 0.0;
        }
        private double dist(int genomeA, int genomeB)
        {
            double marginalOfA = probData.getGenomeMarginalDistribution(genomeA);
            double marginalOfB = probData.getGenomeMarginalDistribution(genomeB);
            double pai1 = marginalOfA / (marginalOfA + marginalOfB);
            double pai2 = 1 - pai1;
            ArrayList conditionalOfGenomeA = probData.getConditionalDistribuationOfGenome(genomeA);
            ArrayList conditionalOfGenomeB = probData.getConditionalDistribuationOfGenome(genomeB);
            ArrayList pGag = calcPGag(conditionalOfGenomeA, pai1, conditionalOfGenomeB, pai2);
            double jsResult = pai1 * (this.Dkl(conditionalOfGenomeA, pGag,1)) + pai2 * (this.Dkl(conditionalOfGenomeB, pGag,2));
            aLocations.Clear();
            bLocations.Clear();
            return (marginalOfA + marginalOfB) * jsResult;
        }

        private ArrayList calcPGag(ArrayList conditionalOfGenomeA, double pai1, ArrayList conditionalOfGenomeB, double pai2)
        {
            int aCount = conditionalOfGenomeA.Count;
            int bCount = conditionalOfGenomeB.Count;
            int i = 0, j = 0;
            ArrayList pGag = new ArrayList();
            while (i < aCount && j < bCount)
            {
                int currentLocation = pGag.Count;
                if ((dCell)conditionalOfGenomeA[i] == (dCell)conditionalOfGenomeB[j])
                {
                    pGag.Add(new dCell(((dCell)conditionalOfGenomeA[i]).CogNo, ((dCell)conditionalOfGenomeA[i]).multBy(pai1) + ((dCell)conditionalOfGenomeB[j]).multBy(pai2)));
                    this.aLocations.Add(currentLocation);
                    this.bLocations.Add(currentLocation);
                    i++;
                    j++;
                }
                else if ((dCell)conditionalOfGenomeA[i] > (dCell)conditionalOfGenomeB[j])
                {
                    pGag.Add(new dCell(((dCell)conditionalOfGenomeB[j]).CogNo, ((dCell)conditionalOfGenomeB[j]).multBy(pai2)));
                    this.bLocations.Add(currentLocation);
                    j++;
                }
                else if ((dCell)conditionalOfGenomeA[i] < (dCell)conditionalOfGenomeB[j])
                {
                    pGag.Add(new dCell(((dCell)conditionalOfGenomeA[i]).CogNo, ((dCell)conditionalOfGenomeA[i]).multBy(pai1)));
                    this.aLocations.Add(currentLocation);
                    i++;
                }
            }
            if (j < bCount)
            {
                while (j < bCount)
                {
                    this.bLocations.Add(pGag.Count);
                    pGag.Add(new dCell(((dCell)conditionalOfGenomeB[j]).CogNo, ((dCell)conditionalOfGenomeB[j]).multBy(pai2)));                 
                    j++;
                }
            }
            if (i < aCount)
            {
                while (i < aCount)
                {
                    this.aLocations.Add(pGag.Count);
                    pGag.Add(new dCell(((dCell)conditionalOfGenomeA[i]).CogNo, ((dCell)conditionalOfGenomeA[i]).multBy(pai1)));
                    i++;
                }
            }
            return pGag;
        }

        private double Dkl(ArrayList p, ArrayList pGag,int t)
        {
            double dklResult = 0;
            int pCount = p.Count;
            if (t==1)
                for (int i = 0; i < pCount; i++)
                    dklResult = dklResult + ((dCell)p[i]).Length * Math.Log(((dCell)p[i]).Length / ((dCell)pGag[(int)aLocations[i]]).Length, 2.0);
            else
                for (int i = 0; i < pCount; i++)
                    dklResult = dklResult + ((dCell)p[i]).Length * Math.Log(((dCell)p[i]).Length / ((dCell)pGag[(int)bLocations[i]]).Length, 2.0);

            return dklResult;
        }
    }
}